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ABSTRACT 

A  bulk  turbulence-closure  mixed  layer  model  is  generalized  to  allow  prediction  of  very 
deep  polar  sea  mixing.  The  model  includes  unsteady  three-component  turbulent  kinetic 
energy  budgets.  In  addition  to  terms  for  shear  production,  pressure  redistribution,  and 
dissipation,  special  attention  is  devoted  to  realistic  treatment  of  thermobaric  enhancement 
of  buoyancy  flux  and  to  Coriolis  effects  on  turbulence.  The  model  is  initialized  and 
verified  with  CTD  data  taken  by  R/V  Valdivia  in  the  Greenland  Sea  during  winter  1993- 
1994.  Model  simulations  show  (i)  mixed  layer  deepening  is  significantly  enhanced  when 
the  thermal  expansion  coefficient's  increase  with  pressure  is  included;  (ii)  entrainment  rate 
is  sensitive  to  the  direction  of  wind  stress  because  of  Coriolis;  and  (iii)  the  predicted 
mixed  layer  depth  evolution  agrees  qualitatively  with  the  observations.  Results 
demonstrate  the  importance  of  water  column  initial  conditions,  accurate  representation  of 
strong  surface  cooling  events,  and  inclusion  of  the  thermobaric  effect  on  buoyancy,  to 
determine  the  depth  of  mixing  and  ultimately  the  heat  and  salt  flux  into  the  deep  ocean. 
Since  coupling  of  the  ocean  to  the  atmosphere  through  deep  mixed  layers  in  polar  regions 
is  fundamental  to  our  climate  system,  it  is  important  that  regional  and  global  models  be 
developed  that  incorporate  realistic  representation  of  this  coupling. 
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I.  INTRODUCTION 
A.         BACKGROUND 

All  interaction  of  the  atmosphere  with  the  ocean  must  occur  where  the  two  fluids  are 
in  contact  -  that  is,  in  the  upper,  fully  turbulent  layer  of  the  ocean  known  as  the  mixed  layer. 
This  mixed  layer  is  typically  more  buoyant  than  the  water  below,  so  that  changes  in  the 
atmosphere  are  to  some  extent  isolated  from  the  deep  ocean  by  the  more  buoyant  layer. 
Exceptions  to  this  general  structure  are  in  the  polar  ocean  regions,  where  mixing  may  occur 
to  great  depths,  or  even  to  the  bottom.  This  coupling  of  the  ocean  to  the  atmosphere  through 
deep  mixed  layers  in  polar  regions  is  fundamental  to  the  initiation  of  the  global  thermohaline 
conveyor  belt  and  to  the  earth's  climate  system. 

Mixed  layer  modeling  has  historically  used  a  linear  dependence  of  density  upon 
temperature,  so  that  the  thermal  expansion  coefficient  is  treated  as  a  constant.  This 
approximation  is  valid  over  much  of  the  ocean,  where  mixed  layers  are  typically  shallow,  but 
breaks  down  in  the  ocean  regions  where  there  are  large  excursions  in  mixed-layer  depth 
and/or  mixed-layer  temperature.  The  equation  of  state  for  seawater  shows  that  the 
dependence  of  density  and  buoyancy  upon  temperature  is  not  simply  linear;  the  thermal 
expansion  coefficient  varies  with  pressure,  and  does  so  most  at  temperatures  near  freezing 
(Garwood,  1991).  The  thermobaric  effect  on  the  stability  of  a  water  column  is  most 
significant  then  where  nearby  parcels  -  very  cold  and  at  depth  -  have  strong  contrast  in  their 
potential  temperatures.  This  is  a  characteristic  state  common  to  deep  mixed  layers  in  the 
polar  seas,  where  traditional  mixed  layer  modeling  has  failed  to  predict  realistic  mixing  and 
deep  water  formation. 
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Another  common  approximation  that  may  not  be  valid  for  deep  mixed  layers  is  the 
neglect  of  planetary  rotation's  effects  on  turbulence.  In  more  shallow  mixed  layers,  this 
neglect  is  reasonable  because  the  scale  of  the  turbulence  is  much  smaller  than  the  distance 
over  which  the  effects  of  rotation  can  be  felt,  the  Rossby  radius.  However,  in  polar  regions 
the  Rossby  radius  is  much  reduced  (~1  km).  The  scale  of  turbulent  eddies  is  equivalent  to 
the  depth  of  mixing,  so  if  deep  mixing  occurs  (~1  km),  the  eddies  are  large  enough  to  be 
affected  by  planetary  rotation,  and  the  effects  need  to  be  included  in  models  of  deep 
entrainment.  Furthermore,  the  vector  planetary  rotation  (q  =  q  e  +q.§z)  needs  to  be 
considered  because  of  the  three-dimensionality  of  turbulence. 
B.         OBJECTIVE 

Deep  mixing  is  controlled  by  factors  more  complex  than  those  of  traditional  mixed 
layer  physics;  the  proper  representation  of  non-linear  buoyancy  effects,  energy-converting 
rotational  effects,  non-steady  atmospheric  forcing  and  water  column  initial  conditions  is 
critical  to  calculating  the  fluxes  of  heat,  salinity,  and  tracers  into  the  deep  ocean.  In  this 
research  a  model  was  developed  and  tested  that  expands  on  traditional  mixed  layer  modeling 
by  adding  the  effects  of  thermobaricity  and  planetary  rotation,  that  allows  initialization  with 
hydrographic  data,  and  that  forces  the  mixing  with  realistic  meteorological  time  series  of 
surface  heat  flux  and  wind  stress.  The  goal  is  not  just  to  compare  the  final  mixed  layer 
depth  to  that  observed,  but  to  see  how  the  observed  event  is  simulated  in  terms  of  relative 
contributions  of  energy  and  evolution  in  time. 

Chapter  II  summarizes  the  state  of  the  art  leading  up  to  this  research.  The  model 
equations  are  developed  in  Chapter  III,  and  the  numerical  approach  is  discussed.  Chapter  IV 


presents  the  oceanographic  data  used  to  initialize  and  the  meteorological  data  used  to  force 
the  model.  Chapter  V  shows  model  results  for  hypothetical  cases  with  constant  forcing, 
followed  by  results  with  time  dependent  forcing.  Chapter  VI  summarizes  the  results  of  the 
research  and  makes  conclusions  and  recommendations  for  further  study. 


II.  MIXED  LAYER  THEORY 

A.    GENERAL  MIXED  LAYER  AND  DEEP  WATER  FORMATION  THEORY 

Mixed  layer  theory  concerning  normal  mixing  depths  (<100  m)  has  been  well 
developed  and  is  summarized  by  Zilitinkevich  et  al.  (1979),  Garwood  (1979),  and  more 
recently  by  Large  et  al.  (1994).  However,  none  of  these  has  direct  application  to  very  deep 
mixing  (order  1  km  and  greater). 

Early  speculations  that  deep  water  is  formed  not  just  by  drainage  off  continental 
shelves,  but  by  convection  in  the  open  ocean,  were  made  by  Nansen  (1 912)  for  the  Arctic  and 
by  Wiist  (1928)  for  the  Antarctic;  it  was  not  until  much  later  that  there  was  enough  evidence 
to  prove  that  open  ocean  deep  convection  does  occur. 

Brown  and  Beardsley  (1978)  suggested  that  overturning  during  winter  storms  could 
be  the  primary  source  of  mixing  that  forms  Gulf  of  Maine  Intermediate  Water.  They  were 
cautious  at  the  time,  however,  because  they  felt  their  data  did  not  completely  prove  this 
result. 

Gordon  (1978)  reported  the  discovery  of  a  very  localized  region  (-13  km)  of 
convection  in  the  Weddell  Sea.  Killworth  (1979)  showed  that  the  entire  region  of  the 
Weddell  Sea  gyre  was  susceptible  to  overturning,  and  that  reduced  stratification  within 
cyclonic  eddies  could  be  the  mechanism  that  "selects"  such  narrow  areas  of  convection 
within  the  larger  gyre. 


B.         THERMOBARICITY 

Gill  (1973),  considering  the  mixing  of  water  flowing  off  the  continental  shelf  into  the 
Weddell  Sea,  noted  the  importance  of  the  variation  of  the  thermal  expansion  coefficient  with 
depth  in  cold  water.  He  included  this  effect  in  his  calculations  and  determined  that  there 
could  be  parcels  of  water  with  potential  temperature  -  salinity  (  e  -s)  properties  such  that, 
if  displaced  a  finite  distance  downward,  the  parcels  would  be  unstable,  their  relative  density 
to  the  surrounding  water  continuing  to  increase  as  they  descended,  thus  accelerating  the 
parcels  downward. 

McDougall  (1987)  coined  the  term  thermobaricity  for  the  increase  in  seawater's 
thermal  expansion  coefficient  with  pressure  (p  ),  or  the  —  (  — £  )  term  in  the  equation  of 
state,  investigated  its  effect  on  mixing  across  neutral  surfaces,  and  discovered  that  although 
both  upward  and  downward  vertical  velocities  across  the  neutral  surface  are  possible,  they 
are  usually  downward. 

Garwood  (1991)  studied  thermobaricity' s  effect  on  mixing  across  a  non-neutral 
surface,  the  interface  at  the  bottom  of  the  mixed  layer.  Although  previously  neglected,  the 
buoyancy-flux  enhancement  by  thermobaricity  was  found  to  be  an  important  source  of 
vertical  turbulent  kinetic  energy,  leading  to  increased  entrainment. 

Garwood  et  al.  (1994)  noted  that  in  addition  to  enhancing  entrainment  rates, 
thermobaricity  might  assist  in  bringing  surface  water  to  depth  by  causing  water  column 
instabilities.  A  parcel-type  instability  would  occur  when  a  parcel  from  the  cold  and  fresh 
turbulent  layer  overshoots  the  interface,  and  continues  to  accelerate  downward  because 
pressure  compresses  it  more  than  the  surrounding  water  is  compressed.    A  layer-type 


instability  would  occur  as  an  interface  advected  downward  becomes  unstable  due  to  the 

thermobaric  effect. 

C.         ROTATION  IN  ATMOSPHERIC  AND  OCEANIC  BOUNDARY  LAYERS 

Rossby  and  Montgomery  (1935)  were  the  first  to  suggest  that  geophysical  boundary 
layers  should  scale  vertically  with  — -  ,  where  u,  is  the  friction  velocity  and  f  =  2q3  is  the 
vertical  component  of  planetary  rotation.    Hence  a  boundary  layer  "Rossby  number"  is 

~h~f~ 

Gascard  (1973)  showed  that  in  a  homogeneous  water  column,  large  amplitude 

internal  waves  generated  by  wind  surges  in  shallow  regions  can  propagate  upward  because 

of  planetary  rotation,  bringing  energy  with  them,  and  suggested  they  may  couple  with  a  and 

be  involved  in  mixing  and  deep  water  formation. 

Wyngaard  et  al.  (1974  )  were  the  first  to  evaluate  turbulence  terms  containing  q2  ; 

their  study  was  of  the  convective  atmospheric  boundary  layer,  and  found  q2  's  effect  to  be 

small.    However,  turning  to  ocean  turbulence,  Garwood  (1977)  found  a  Rossby  number 

dependence  for  oceanic  mixed  layers  in  a  bulk  turbulence  closure  model.   Garwood  et  al. 

(1985)  went  on  to  show  that  planetary  rotation  acts  to  exchange  turbulent  kinetic  energy 

between  horizontal  and  vertical  components  in  a  steady-state  mixed  layer.  Vertical  turbulent 

energy  and  thus  mixed  layer  depth  are  increased  under  easterly  winds.   Garwood  (1991) 

extended  this  hypothesis  to  the  more  general  case  of  deepening  mixed  layers. 


III.  MODEL 

A.         DERIVING  THE  EQUATIONS  FOR  TURBULENT  KINETIC  ENERGY 

Representing  the  Navier-Stokes  equations  of  motion  as 

1    dp 


d  u .  d  u ..  i     3  n  d2  u . 


u 


2e..,Q.u,     -  5.,g    +v i  (1) 


at  j  axj       p  ax.  «*  j  *       i3  dxi 

the  quantity  turbulent  kinetic  energy  is  formed  by  decomposing  the  variables  into  their  mean 
and  turbulent  parts,  for  example,  representing  one  of  the  three  velocity  components  as 

"i  =  ~i  +  ui  (2) 

where  the  turbulent  part  is  not  assumed  to  be  small  compared  to  the  mean.  Here  i  =  1,2,3 
corresponds  to  the  Cartesian  coordinate  directions  of  east,  north,  and  up; 
[  q_.  ]  =  q  [  0  cos  0  sincf  ] ;  molecular  viscosity  is  v ;  gravitational  acceleration  is  g;  and 
<p  is  latitude.  The  u-component  kinetic  energy  per  unit  mass  contained  in  the  turbulence  will 
then  be 


u     TKE  =   -^—  (3) 

1  2 


Since 


g         u{  d  u,        d  u 

at  (~2_)  =  Ll  at 


u,  — —   -   u 


(4) 


and  the  right-hand-side  contains  only  quantities  predicted  by  the  equations  of  motion,  kinetic 
energy  can  be  calculated.  Again  generalizing  to  the  three  components,  the  last  term  can  be 


written  as 
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after  Reynolds'  decomposing,  averaging,  and  making  the  Boussinesq  approximation  on 
equation  (1).  The  first  term  on  the  right-hand-side  of  equation  (4)  is  formed  by  performing 
the  multiplication  and  decomposition  first,  and  then  averaging,  resulting  in 


du . du  du±  u^du±  du. du 
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Subtracting,  we  have 
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ax; 
The  first  and  third  terms  on  the  right-hand-side  are  simplified  by  using  the  reverse  product 

rule  for  differentiation  and  the  incompressible  form  of  the  continuity  equation.  Then 


dul  id 
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and 
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The  TKE  equation  for  the  system  is 
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Changing  notation  so  that  [  u.]  =  [  u,  v,  w  ] ,  [xi]  =  [  x,  y,  z  ]  and  assuming 

horizontal  homogeneity  so  that  — -  (       )   =  — —  (       )   =     gives  the  component  equations 

ox  ay 

for  turbulent  kinetic  energy: 


3u'  „ du  d     .  2   v  2         3u 

2u'w——     -   — —  (   w  u  ""   )      +   — pr 


3  t  3z  3z  pn         5x 


4Q2v  w     +  4Q3u'  v 


,         32  u '              d2  u '              d2  u '  / 1 1  \ 

v  (   v +  w +   w (11) 

dx2  dy2  dz2 


3v2  _ dv  d     .  2   v  2      ,  dv  ._  ■ 

- —     =     -  2v  W  — -     -   -—  (   w  v2   )      +  — p'  — —     -  4Q 
dt  dz  dz  p  dy 


d2  v '  d2  v '  d2  v 

+  V     (    V"  +    V  +    V 


dx2  dy2  dz- 
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(12) 


dw2  d     .  —3-  1      d  — „    1         dw  .    g 

— —     =     -   ■=—  (   w     )      -  2 —  w p'     +2  — p'  — —     -  2-2-p'  w 

ot  dz  p     dz  p  d  2  p. 


5  w"  d2  w  d2  w 


+   4Q0u'W      +v    (    w '  —    +   w '  —    +   w  •  )  ( 1  J) 

dx2  dy2  dz2 

where  the  terms  on  the  right-hand-side  are  shear  production,  turbulent  transport,  pressure 

redistribution,  pressure  transport,  buoyant  production/damping,  Coriolis  redistribution,  x-z 

rotational  redistribution,  and  viscous  dissipation;  and  the  term  on  the  left-hand-side  is  called 

the  storage  term. 

B.         INTEGRATING  FOR  BULK  TKE 

Integrating  over  the  depth  of  the  mixed  layer  plus  entrainment  zone  will  give  the 

bulk  TKE  for  the  layer.  The  storage  term  becomes 

f°  ?Hdz  =   dh<~2>  (14) 

J  -h      d  t  dt 

where  <  >  denotes  the  vertically  integrated  value.  Turbulent  transport  and  pressure 
transport  terms  integrate  simply  to  zero,  assuming  no  turbulent  energy  or  pressure 
perturbations  are  transported  across  the  boundaries.  (Entrainment  does  not  cross  the 
interface  at  z  =  -h;  rather,  it  changes  the  depth  of  h.) 


dw  u" 


f°    -^-^dz  =  wu-2\Q  -ww2\    .     =  0  (15) 

J  -n  a  z 


f°  -A^_Pldz  =  „.p.|    -w-p<U    =0  (16) 

J-h       p0       dz  ° 

Integration  of  the  buoyant  production/damping  term  and  the  x-z  rotational  redistribution  term 
requires  consideration  of  the  approximately-linear  flux  profiles  of  temperature,  salinity,  and 
momentum  within  the  layer,  and  the  "jump  condition"  at  the  bottom  of  the  layer.  For  any 
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conserved  quantity  C,  the  linear  flux  profile  has  a  shape  such  as  Figure  (1)  and  the  jump 
condition  is 


C'W  \_b   =   -weAC  (17) 

This  introduces  a  new  variable  we,  the  entrainment  velocity.  In  the  absence  of  upwelling  or 
downwelling, 


dh 

— —   -  w 

dt 


Alternatively,  if  non-zero  mean  vertical  velocity  is  included,  then 


dh 
dt 


(18) 


(19) 


The  temperature  flux  at  the  surface  is 


Q, 


T'W\Q  =  -^-  (20) 

where  q0  is  the  rate  of  cooling  of  the  surface,  and  c  is  the  specific  heat  of  seawater  at 
constant  pressure.  The  surface  salinity  flux  is 


S'w|0  =  s  vol  (21) 

where  Vol  is  the  rate  of  addition  or  removal  of  water  from  the  surface,  the  net  result  of 
precipitation,  evaporation,  melting,  and  freezing,  expressed  in  units  of  length  per  time.  The 
equations  for  the  fluxes  as  functions  of  z  then  become 

¥•  w  (z)    =   L_(l  +  -£)    -WpAT(-^)  (22) 

and 


Pn  c  h  h 


S'w(z)    =  S  Vol  {.!  +  —)    ~wAS( )  (23) 

h  e  h 

Rather  than  treating  the  thermal  expansion  coefficient  as  a  constant,  it  is  allowed  to  vary 
linearly  with  depth. 

a  (z)    -  a0-  axz  (24) 
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►     WW 


WW 


Figure  1 .  Flux  profiles  shown  are  for  momentum,  but  would  have  similar  shapes  (piecewise 
linear)  for  any  conserved  quantity.  The  upper  plot  is  for  wind  direction  west  to  east,  and  the 
lower  plot  is  for  wind  east  to  west. 
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Then  the  integration  of  the  buoyant  production  and  x-z  rotational  redistribution  terms  reduces 
to  evaluating  definite  integrals  of  polynomials,  so 


f°   -2-2-p'  w   dz  =2  f°  b'  w   dz 

J -h  D„  J  -h 


=      2 


f°  ol(z)    T'W   dz  -  2  f°   3  S'W   dz 

■I  -h  J  -h 


QQ  ah 

gh  [   — °-  (   a.  +   — L-   )    +  3  S  Vol 
pc  3 


a,h 
1.5 


ghwj  AT(  a0+  ^   )    -  3  AS  ]  (25) 


and 


[°    4Q,u'Wdz     =     2Q,h(    ^    +   wAC/)  (26) 

By  scaling,  the  second  term  on  the  right-hand-side  of  equation  (26)  is  small  and  is  neglected 
hereafter. 

The  integral  of  the  shear  production  term  has  contributions  only  from  the  regions 
where  there  is  velocity  shear,  near  the  surface  and  just  below  the  bottom  of  the  fully  turbulent 
layer,  as  shown  in  Figure  (2).  Measurements  have  shown  the  velocity  shear  just  below  the 
surface  to  be  proportional  to  the  friction  velocity,  so  that 

6u  =  iruu.sincp  (27) 


where  <p  is  the  angle  from  north  clockwise  to  the  wind  vector,  u2,  =  — Jx2x  +  \\  ,  and  m. 
is  the  constant  of  proportionality.  The  momentum  fluxes  at  the  surface  and  within  the 
entrainment  zone  are  approximated  by  the  piecewise  linear  functions  in  Figure  (1). 
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Figure  2.  To  model  the  production  of  turbulence  by  shear  in  the  mean  flow,  values  of 
vertical  shear  at  the  ocean's  surface  and  in  the  entrainment  zone  are  estimated  rather  than 
assumed  infinite. 
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The  surface  momentum  flux    ww  1 0  and  the  shear  near  the  surface  —  1 0  are  always 
opposite  in  sign,  so  that  their  product  is  always  negative.  Using  the  relationship 

—  =  u;  sin$  (28) 

the  integral  becomes 


2  1  u'W    — —  dz     =   -2  1       ww    -— dz     -2  1         ww    —— dz 

J-h-b  OZ  J-bz  OZ  J-h-b  OZ 


2   ww  |.  bu     -2   ww  I        s   Au 

u  -/]-  — 

2 


i^i(iii)2  +  ,(Ai)  (29) 


Considering  the  Coriolis  redistribution,  the  Reynolds'  stress  is  parameterized  as 
proportional  to  the  horizontal  shear  of  the  mean  flow,  which  under  horizontal  homogeneity 
is  zero.  This  means  that  because  both  forcing  and  response  are  the  same  everywhere,  there 
can  be  no  net  transport  of  u  turbulence  in  the  v  direction,  and  no  net  transport  of  v  turbulence 
in  the  u  direction. 

4Qj°  w~^dz    «  (  |I  +  |Z  )     =o  (30) 

i J  -h  oy        ox 

Using  the  "return  to  isotropy"  assumption  of  Rotta  (1951),  the  pressure  redistribution 
term  becomes 


3 


J-^P'lRLdz     =  2m7{   <E>  -  3<u-2>   )<E>2  (31) 

P0J-h  dx 

where  <£>  =  (  <~w2~>  +  <~>  +  <w2>  )  and  m2  is  a  constant  of  proportionality. 
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Applying  the  product  and  chain  rules  for  differentiation,  the  dissipation  integral  has 
terms  of  two  forms, 


r o  ,        6  u'  a"  w 

+  V  /  U  '  +    U  '  +    u 

3-h  dx2  dy2 


dy2  dz- 


ro    ,    dv    dv         dv    dv         dv    dv    ,    , 

v  /       ( + + )  dz 

J -h         dx     dx  dy     dy  dz     dz 


f°V2(^_)dZ  (32) 


The  last  term,  viscous  diffusion  of  turbulent  energy,  is  negligible  compared  to  both  buoyant 
and  shear  production,  or  to  turbulent  transport.  The  remaining  terms  on  the  right-hand-side 
of  equation  (32)  are  u-component  viscous  diffusion  and  represent  all  of  the  dissipation  of  u 
TKE.  Total  dissipation  is  assumed  to  be  a  function  of  the  net  turbulent  energy  only,  and  is 
thus: 

D  =  2m,<E>2  (33) 

where  m1  is  a  constant  of  proportionality  and  <  >  denotes  the  mixed-layer  average  value. 
Since  viscous  dissipation  occurs  at  the  smallest  scales,  where  kinetic  energy  is  converted  to 
heat,  total  dissipation  is  also  assumed  to  be  isotropic.  The  corresponding  terms  of  each 
component  equation  would  then  be  one  third  of  the  total  dissipation. 

Using  the  integral  results  of  equations  (14),  (25),  (26),  (29),  (31),  and  (33),  the  bulk 
equations  for  the  components  of  turbulent  kinetic  energy  are: 
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dh<  n ,2  >  2m,       t       ,  t 

3   (_JS)2    +    w    (Au)2    +   2£2,h-S 


at  u.     pc 


2  ■ 


3 
.2.  ,."  2 


+  2ra2(   <£>   -  3<u'2>   )  <E>  2    -   ±m1<E>2  (34) 


dh<v2>  2m        i 

=- =     (  -*■ )      +   w    (  A  v , 

at  u  p„ 


+  2  m,  (  <£>  -  3<v2>  )  <£>  2  -  —  m  <E>2  (35) 


at  pc    u    3 

p 


ghw  [  AT(  a    +  — i—  )  -  3  AS  ]   -  2Q,h  — 
e       °   1 . 5  2   p„ 


+  2m2  (   <£>  -   3<w'2>   )  <E>2    -   —ml<E>2  (36) 


C.         TEMPERATURE  AND  SALINITY  EQUATIONS 

Conservation  principles,  plus  horizontal  homogeneity,  allow  the  potential  temperature 
and  salinity  equations  to  be  written  as 


ae  -36  99'  w                                                                                            rms 

T7  =     "  w^~-     ~  — 5 C37) 

at  oz  dz 

and 


5S               —  3S          3S' v  /"5o\ 

=     -  w  -         (J8) 

at  az  az 

Integrating  over  the  mixed  layer  plus  entrainment  zone,  and  applying  the  jump  condition, 

equation  (17),  the  equations  become 
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|i  -  1  <  -°S-     -  .eA9  )                                                                           (39) 

at  h     p0  cp 

and 

il  =   1  (    SVol     -   w  AS   )                                                                                                     (40) 

at  h                                                                                    y    ' 


D.         BULK  MOMENTUM  EQUATIONS 

Equation  (1),  after  decomposing,  averaging,  making  the  Boussinesq  approximation, 
assuming  horizontal  homogeneity,  and  neglecting  viscous  effects  on  the  mean  flow,  gives 
the  component  equations  for  mean  flow: 


du               —  du            1    dp         ~~  —         -^  —         dvw  /-ai\ 

— -     =     -  w—     -   — -£      -2Q,w     +  2Q  v     -   — (41) 

at  dz  p0   dx  3  dz 


dv  —  dv  1    dp  or^~  dvw 


-  w 


_^i^     -  2Q,u     -   ^-^-  (42) 


dt  dz  pQ   dy  dz 


dw  —  dw  1    dv  ,   p    po  .  «_  —  3ww  (ai\ 

a-'  =     "  w-—     -  — -£-     -g( )      +  2Q,w     -   — (43) 

at  oz  p0   dz  pc  dz 

The  hydrostatic  approximation  is  made,  so  that  equation  (43)  becomes 

|£    =  -Pg  (44) 

oz 
and  upwelling  is  neglected  here. 

Although  the     q2  terms  were  retained  in  the  turbulence  equations  because 

w  ~  o(w) ,  they  can  be  dropped  here  because   w  «  0(11,17)  .  Also  if  the  geostrophic  and 

wind-driven  flows  are  separated,  equations  (41)  and  (42)  become 


dU              —  dU~^—         dvw  (az\ 

—     =     -   w—     +  2Q   V    -   — (45) 

dt  dz  i  d  z 


and 
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dv  -dV  -         dvw 

—     =     -   w—      -  2Q    [/     -   — — 
at  dz  dz 


2Q3T/    -  ^LJ^-  (46) 


where  u  and  v  represent  the  wind-driven  flow  averaged  over  some  short  period  of  time. 
Integrating  equations  (45)  and  (46)  yields 


*(*<"»  =2Q3h<V>    +1*  (47) 

at  3  Po 

and 

*tna<v>j  =  -  2Q3h<a>    +  i  (48) 

3t  3  pn 


E.         ENTRAINMENT  VELOCITY 

Entrainment  velocity  must  now  be  parameterized  to  close  the  system  of  equations; 
the  local  TKE  budget  at  the  bottom  of  the  mixed  layer,  where  the  entrainment  is  occurring, 
can  give  a  scale  estimate.  The  sum  of  equations  (1 1),  (12),  and  (13)  retains  contributions  to 
the  storage  term  from  shear  production,  turbulent  transport,  buoyant  damping,  and 
dissipation: 

du       dv   , 


-      =  -   2  (    V  W  -r—    +    V  W  —— 

at  a  z  a z 


az  p0 


2-2-p-w     -e  (49) 
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i_                                                d  u '              d  u '  d~  u ' 

where         e  =      +  v  (  u  •  —      +  u  • +  u  ■ 


dx2  dy2  d z2 


dv              dv'  d~  v ' 

v'  +  v" +  v 


dx2  dy2  dz: 


dW  d~W  d~  W      ,  srr\\ 

w +   w —    +   w —    )  (-)^) 

dx2  dy2  d z2 

Near  z  =  -h  both  shear  production  and  dissipation  are  small  compared  to  the  turbulence 
transported  from  above  and  to  the  amount  being  removed  by  buoyant  flux,  so  that  the  balance 
is  largely  among  storage,  turbulent  transport  and  buoyant  damping. 


d  (    u'~  +  v"-   +  w'   )    _  d  p- 


=   -   -^-  [   w  (E  +  2^^)    ]-  2-2-p'W  (51) 


d  t  dz  pn  p 


Each  of  these  terms  can  be  scaled  on  bulk  variables,  and  the  time  scale  for  changes  in 
turbulent  energy  at  the  bottom  of  the  mixed  layer  can  be  thought  of  as  the  time  required  for 
a  change  in  energy  content  to  be  turbulently  diffused  from  the  surface,  which  is  posed  here 


to  depend  upon  h  and  <w2> . 

d{  IP2  +  ~  +  ^r2   )  <E>we 


dt 


(52) 


^_[  „•  (B+  2^_)   ]  -  -J 1 (53) 

dz  p0  h 


-2-2-p-w     =     -w    gh(aAT-(3AS)  (54) 

Combining  and  solving  for  entrainment  velocity  give 


m    <E>  \J  <w~> 

W     =    — = — ^r- 

<£>  +  gh  (a&T  -  (3AS; 


(55) 
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F.         NUMERICAL  APPROACH 

The  system  of  eight  equations,  summarized  as  Table  (1),  was  cast  as  a  MATLAB 
function,  and  the  function  was  then  called  using  the  MATLAB  utility  ODE45,  a  Runge-Kutta 
ordinary  differential  equation  solver  that  has  automatic  step-size  control.  The  tolerance  for 
the  solver  was  relaxed  from  its  default  value  of  10"6  to  10"4  for  faster  convergence. 
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3B  1     ,           QQ                                yr    x 

—  =     -  (    —  -  v  A6   ) 

dt  h        pnc              e 

as  i    -                - 

=     —  (   SVol  -  w  AS  ) 

at  h                  e 


d  (h<U>)  nn  ^    -  Tx 

— - =    2Q,h<V>     +   — 

dt  3  p 


o 
dt  3  p 


dib<V>)      -    -2Q.h<5>     ,2Z 


0 


w 


m  <E>\j  <w'2> 


<£>   +  gh[  (aQ  +  a^)  AT  -   3  AS] 


dh<u'2>               2m3  ,   Txv2  ,     -v2        «       ^  T 
= (  —  )2   +  w    (Au)2   +  2Q,h- 

3t  «  Pn  P 


°  2 

2^    \^c*-^2  2  -„^2 


+  2m2  (  <£>  -  3<u,4i>  )  <£>"    -   —  m1<E> 


T 


(  —  )      +  wr   (Av) 


at  u      p, 


2<w    \^*.^2     _      2_    ^-.^  2 


+  2^  (  <£>  -  3<v"^>  )<£>'    -   —  m1<E> 


dh<W'2>  Qq  3!h  -    Tr 

— at-    =    gh[  ^(a°+T)+MM 

-  ghwe[  AT(  a0+   -i-   )    -  3  AS  ] 

1 . 5 

-  2Q,h  — 

_         o 

+  2hl(  <E>   -   3<w'2>  )<E>2    -    —  m<E> 
2  3     * 


Table  1 .  Equations  Used  in  the  Model 
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IV.    OBSERVATIONAL  DATA 

A.         IMPORTANCE  OF  THE  GREENLAND  SEA  REGION 

As  summarized  by  Visbeck  et  al.  (1995),  the  Greenland  Sea  has  long  been  considered 
to  have  characteristics  favorable  for  deep  water  formation.  Nansen,  in  1 906,  observed  only 
small  differences  in  properties  of  Greenland  Sea  surface  and  deep  waters,  and  proposed  that 
deep  water  formation  occurred  there.  Similar  (but  more  modern)  observations  were  made 
by  Rudels  et  al.  (1989),  when  mixing  seemed  to  have  recently  extended  to  1500  m.  A 
synthesis  of  the  physical  explanations  provided  by  Visbeck  et  al.  (1995)  and  Schott  et  al. 
(1993)  for  how  the  region  becomes  favorable  to  deep  convection  includes  the  influences  of 
the  previous  seasons  as  well  as  local  forcing.  Bourke  et  al.  (1992)  showed  that  during  the 
ice-free  months  of  the  year,  upper  water  conditions  within  the  basin  can  be  variable  due 
largely  to  advection  by  the  Jan  Mayen  Current.  In  early  winter,  as  the  marginal  ice  zone 
propagates  rapidly  eastward,  brine  release  and  entrainment  act  to  increase  the  salinity  of  the 
surface  waters,  beginning  to  erode  the  buoyancy  of  the  mixed  layer.  Deepening  during  this 
period  may  extend  far  enough  into  the  warm  Atlantic  Intermediate  Water,  so  that  when 
subsequent  cooling  and  stirring  occur,  warming  by  entrainment  will  dominate  over  surface 
cooling;  the  formation  of  ice  is  prevented  and  an  ice-free  area,  surrounded  on  three  sides  by 
ice  but  open  to  the  north,  known  as  the  Nordbukta,  forms.  The  Nordbukta,  situated  over  the 
Greenland  Sea  gyre,  is  then  preconditioned  for  deep  convection.  Hence,  the  data  discussed 
below,  provided  by  Gascard  and  Lherminier  (personal  communication),  were  very  well 
suited  for  use  as  a  test  case  for  the  model. 
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B.         OCEANOGRAPHIC  DATA 

Two  station  soundings  taken  by  the  R/V  Valdivia  at  74.5  °N,  2.5°  W,  during 
February-March  1994  were  examined.  At  this  station,  the  ocean  depth  is  more  than  3600  m. 
For  each  sounding,  the  pressure,  depth,  in  situ  temperature,  potential  temperature,  salinity, 
potential  density  referenced  to  the  surface,  potential  density  referenced  to  2000  m,  sound 
velocity,  and  dynamic  height  were  available.  In  addition  to  the  oceanographic  data,  each 
sounding  had  recorded  the  station  number,  date,  hour,  precise  latitude  and  longitude,  depth 
of  sounding,  sea  surface  temperature,  sea  surface  salinity,  wind  speed  and  direction,  and  air 
temperature.  Figure  (3)  shows  February  16th's  vertical  profiles  of  potential  temperature, 
salinity,  and  potential  density  referenced  to  the  surface.  A  200-m  deep  mixed  layer  of  cold 
fresh  water  overlies  more  saline  and  warmer  water;  the  buoyancy  effects  of  temperature  and 
salinity  are  in  opposition,  but  salinity  dominates.  The  16  February  profiles  show  that  the 
magnitudes  of  the  temperature  and  salinity  gradients  at  the  bottom  of  the  mixed  layer  are 
very  large,  exceeding  0.1  C/m  and  0.01  psu/m,  respectively.  Figure  (4)  shows  profiles  of  the 
same  quantities  for  March  19th,  31  days  later;  no  intermediate  soundings  were  available. 
The  mixed  layer  is  much  deeper  (600  m),  and  the  entrainment  of  warmer  more  saline  water 
has  diminished  the  contrast  between  the  two  water  masses;  the  buoyancy  jump  has  been 
eroded  by  the  mixed  layer  deepening. 

Below  600  m  the  two  soundings  look  very  similar,  especially  if  the  March  sounding 
is  offset  upward  by  about  30-40  m.  This  suggests  a  slow  but  steady  down  welling  speed  of 
about  1  m/day,  which  is  small  compared  to  the  average  entrainment  speed  of  about  13  m/day. 
Comparison  of  heat  and  salinity  content  between  the  two  soundings  gives  an  average  loss  of 
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Figure  3.  Sounding  of  16  February  1994  shows  a  200  m  deep  mixed  layer  overlying  warmer, 
more  saline  water. 
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Figure  4.  Sounding  of  19  March  1994  shows  mixed  layer  has  deepened  to  600  m. 
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130  W/nr  and  no  surface  flux  in  salinity.  The  similarity  of  profiles  at  depth  as  well  as  the 
apparent  conservation  of  salinity  gives  some  indication  that  horizontal  advection  is  small, 
and  that  one-dimensional  modeling  should  be  useful  for  this  case.  As  will  be  seen  in  Chapter 
V,  the  calculated  Eckman  transport  for  the  layer  suggests  one-dimensional  modeling  as  well. 
C.         METEOROLOGICAL  DATA 

Meteorological  data  from  the  European  Center  for  Medium  Range  Weather  Forecasts 
(ECMWF)  were  available  for  February  2nd  through  May  3 1  st.  Parameters  were  six-hourly 
values  of  surface  temperature,  2  m  temperature,  1 0  m  u-component  of  wind,  lOmv- 
component  of  wind,  surface  sensible  heat  flux,  and  surface  latent  heat  flux.  The  values  were 
interpolated  from  nearest-neighbor  2°  grid  points  in  the  ECMWF  analyses.  Figure  (5)  shows 
the  wind  speed  and  associated  wind  stress  for  the  period  February  1 6th  through  March  1 9th, 
and  Figure  (6)  shows  the  surface  latent  and  sensible  heat  fluxes  for  the  same  period.  The 
mean  heat  loss  for  the  period  indicated  by  these  meteorological  estimates  is  1 82  W/m2, 
significantly  larger  than  the  actual  average  heat  loss  measured  in  the  water  column.  This 
difference  may  be  explained  by  lack  of  inclusion  of  net  radiation  or  by  error  attributable  to 
coarse  resolution  in  the  ECMWF  product.  Methods  of  correcting  for  this  difference  are 
discussed  in  Chapter  V.  Three  events  stand  out  as  periods  of  large  surface  forcing,  or  storms; 
the  latest  of  the  three  is  the  largest.  It  would  be  interesting  to  determine  whether  the  mixed 
layer  depth  was  relatively  stable  until  this  final  large  storm  mixed  the  layer  down  to  600  m, 
or  whether  the  deepening  occurred  gradually. 
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Figure  5.   Meteorological  data  from  ECMWF  shows  that  wind  speed  and  wind  stress  are 
highly  variable  over  the  period. 


SSHF 


SLHF 


100 

0 

-100 

CM 

E-200 

| -300 

-400 

-500 


-600 


•V 

"W 

!V| 

I) 

, 

I 

V 

40  60  80 

Julian  Day 


-200 


-250 


40  60  80 

Julian  Day 


SSHF  +  SLHF 


40  60  80 

Julian  Day 


Figure  6.   l  hree  episodes  of  large  surface  heat  loss  stand  out  in  the  sum  of  surface  sensible 
and  latent  heat  fluxes,  from  ECMWF. 
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V.  MODEL  RESULTS 

A.         EFFECT  OF  THERMOBARICITY 

The  general  behavior  of  the  model  was  tested,  comparing  results  in  cases  where 
thermobaricity  is  included  to  cases  where  it  is  neglected.  These  cases  demonstrate  that  the 
increase  in  thermal  expansion  with  pressure  (depth)  does  increase  the  depth  of  mixing. 
Figures  (7)  and  (8)  show  the  deepening  and  the  TKE  components  for  a  25-day  period  with 
constant  forcing  of  southerly  winds  (0.8217  dynes/cm2  wind  stress)  and  moderate  cooling 
(100  W/m2).  Including  the  thermobaric  effect  increased  deepening  by  50%.  Figures  (9)  and 
(10)  show  the  deepening  and  the  TKE  components  for  a  25-day  period  with  the  same  forcing 
except  stronger  cooling  (400  W/m2);  including  the  thermobaric  effect  increased  deepening 
by  33%.  The  point  where  the  vertical  turbulent  energy  exceeds  the  wind-driven  turbulent 
energy  (v-component  in  this  case)  can  be  interpreted  as  a  shift  from  forced-convection 
dominant  to  free-convection  dominant.  The  shift  of  regimes  occurs  earlier  when 
thermobaricity  is  included. 

Parcels  above  and  below  a  typical  Arctic  region  mixed-layer  interface  have  a  smaller 
buoyancy  contrast  with  al  *  0  than  they  have  with  thermobaricity  neglected;  they  also  have 
a  smaller  buoyancy  contrast  at  depth  than  they  would  at  shallower  mixed  layer  depths.  This 
diminished  buoyancy  jump  will  act  to  increase  we,  increasing  the  deepening  rate.  Further, 
the  buoyancy  flux  at  the  bottom  of  the  layer  is  proportional  to  both  the  entrainment  velocity 
and  the  buoyancy  jump;  entrainment  buoyancy  flux  is  increased  by  both  of  these  factors. 
Another  process  that  deepens  more  due  to  thermobaricity  is  surface  cooling;  heat  taken 
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Figure  7.  With  100  W/m2  of  cooling,  the  model  deepens  the  mixed  layer  50%  more  when 
thermobaricity  is  included  in  the  calculation  (a,  *  0).  Wind  stresses  Taux  and  Tauy  are  in 
dynes/cm2;  heating  QO  is  in  W/m2;  and  a,  is  in  (cm  °C)"'. 
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Figure  8.  Velocity  variances  (  u  -1 ,  for  example)  represent  the  turbulent  energy  present  in 
the  eddy  field  of  each  component.  Comparison  of  the  components  in  the  plots  above  shows 
that  the  increase  in  deepening  for  the  case  with  thermobaricity  is  due  to  increased  turbulent 
energy  in  the  vertical  component. 
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Figure  9.  With  400  W/m2  of  cooling,  the  model  deepens  the  mixed  layer  33%  more  when 
thermobaricity  is  included  in  the  calculation  (a1  *  0).  Wind  stresses  Taux  and  Tauy  are  in 
dynes/cm2;  heating  QO  is  in  W/m2;  and  a,  is  in  (cm  °C)"'. 
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Figure  10.  Turbulent  energy  components,  represented  by  velocity  variances,  show  that  the 
increase  in  deepening  is  due  to  more  vertical  turbulent  energy.  The  shift  from  forced 
convection  to  free  convection  occurs  earlier  with  thermobaricity  included. 


out  at  the  surface  will  have  a  greater  density  effect  at  depth  than  at  shallower  depths. 
Thermobaricity  increased  the  deepening  rate  even  for  cases  with  as  little  as  1 0  W/m2  of 
cooling. 

B.  EFFECT  OF  COOLING 

To  examine  the  sensitivity  to  seasonal  cooling,  the  model  was  integrated  for  a  period 
of  1 00  days  from  the  initial  profile  of  1 6  February,  using  forcing  corresponding  to  the  means 
of  u-wind  and  v-wind  and  two  different  rates  of  cooling.  In  the  first  case,  the  observed 
average  heat  loss  of  130  W/m2  was  used;  in  the  second  case,  twice  that.  Figure  (11)  shows 
that  doubling  the  rate  of  cooling  has  the  effect  of  doubling  the  deepening  rate;  all  other  initial 
conditions  and  forcing  are  the  same,  and  the  heat  is  simply  removed  more  quickly.  Figure 
(12)  shows  that  free  convection  dominates  much  sooner  with  the  stronger  cooling.  It  is  also 
interesting  to  see  the  response  of  the  horizontal  turbulent  energy  components  to  the  dramatic 
increase  of  vertical  turbulent  energy.  They  increase  as  a  result  of  pressure  redistribution;  the 
turbulence  is  attempting  to  become  more  isotropic. 

An  interesting  question  is  whether  it  was  possible  for  the  water  column  as  measured 
on  16  February  to  be  mixed  to  the  bottom  before  summer  ice  melt  and  warming  shut  off  deep 
mixing.  Looking  at  Figure  (11)  again  with  this  question  in  mind,  it  seems  that  there  could 
not  be  enough  time  left  even  with  very  strong  cooling  for  mixing  to  the  bottom  to  have  taken 
place  in  1994. 

C.  EFFECT  OF  Q2  ROTATION 

The  direction  of  the  wind  has  a  potential  effect  on  deepening  because  it  can  oppose 
or  coincide  with  the  direction  of  planetary  rotation.  To  examine  this  effect,  two  1 00-day 
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Figure  1 1.  Doubling  the  cooling  rate  doubles  the  rate  of  deepening. 
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Figure  12.    The  vertical  component  of  turbulent  energy  is  dominant  much  sooner  when 
cooling  is  260  W7m:  vice  130  W/m2. 
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cases  were  simulated,  one  with  large  westerly  wind  stress  (+2  dynes/cm2)  and  one  with  large 
easterly  wind  stress  (-2  dynes/cm2).  Figure  (13)  shows  that  the  westerly  case  deepened  to 
1289  m,  while  the  easterly  case  deepened  further,  to  1468  m.  In  Figure  (14),  while  the  total 
turbulent  energy  for  the  two  cases  is  about  the  same,  the  amount  of  energy  that  is  in  the 
vertical  TKE  is  strikingly  different.  With  easterly  wind,  u  turbulent  energy  is  converted  to 
w  turbulent  energy  (Figure  (15) );  with  westerly  wind,  w  turbulent  energy  is  diminished  by 
conversion  to  u  turbulent  energy  (Figure  (16) ).  The  lower  total  energy  for  the  deeper  case 
is  explained  by  the  fact  that  there  has  been  more  buoyant  damping  (sometimes  called  buoyant 
consumption)  associated  with  the  increased  entrainment. 

D.  CONSTANT  AVERAGE  FORCING  RESULTS  COMPARED  TO  DATA 

When  the  model  is  integrated  for  3 1  days  with  the  mean  wind  components  and 
observed  average  cooling,  a  mixed  layer  depth  is  calculated  that  is  very  close  to  that  observed 
(627  m).  However,  averaging  the  vector  wind  components  under-represents  the  amount  of 
stirring;  in  order  to  make  a  more  valid  comparison,  mean  wind  speed  was  computed  and  a 
typical  wind  direction  chosen.  Components  of  wind  stress  were  then  calculated  from  these, 
and  the  depth  of  mixing  after  integrating  for  31  days  was  733  m  (Figure  (17) ). 

E.  TIME  DEPENDENT  FORCING 

The  model  was  solved  for  the  3 1  days  between  the  times  of  the  CTD  casts  of  1 6 
February  and  1 9  March  by  linearly  interpolating  the  winds  and  heat  fluxes  for  times  between 
the  six  hourly  estimates  from  ECMWF.  A  final  depth  of  802  m  was  reached,  and  Figure  (18) 
shows  that  the  periods  of  more  rapid  deepening  correspond  to  the  more  intense  forcing 
events.  The  deepening  of  the  mixed  layer  did  not  occur  in  one  large  surge  during  the  third 
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Figure  13.  Westerly  winds  (from  the  west,  upper  plot)  are  not  able  to  deepen  the  mixed  layer 
as  much  as  easterly  winds  (lower  plot).  Wind  stress  is  in  dynes/cm2;  heating  is  in  W/m2. 
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Figure  14.  Total  turbulent  energy  is  nearly  the  same  for  the  two  cases,  but  the  amount  that 
is  in  the  vertical  eddies  is  different.  Wind  from  the  east  (lower  plot)  puts  more  energy  into 
the  vertical  component. 
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Figure  15.  With  forcing  by  easterly  winds,  the  contribution  from  rotation  is  negative  for  u 
TKE  and  positive  for  w  TKE:  w  receives  energy  from  u. 


44 


u  TKE  Budget 


50 


< 

w 

n       0 

E 
o 


-50 


2  4  6 

seconds 


x  10 


10 

6 


-  SP-sfc 
-  SP-intfc 

-  p  redistr 
rotation 

-  dissipation 


50 


w  TKE  Budget 


CO 

< 

w 

m   0 

E 
o 


-50 


-  '—^  '  -'  '   '  ■*-■■■.' 


4        6 
seconds 


x  10 


10 

6 


SP-sfc 
SP-intfc 
p  redistr 
rotation 
dissipation 


Figure  1 6.  With  forcing  by  westerly  winds,  the  contribution  from  rotation  is  positive  for  u 
TKE  and  negative  for  w  TKE:  w  gives  up  energy  to  u. 
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Figure  1 7.  The  upper  plot  shows  deepening  over  3 1  days  with  forcing  computed  from  mean 
wind  components.  The  lower  plot  shows  deepening  over  3 1  days  with  forcing  computed  from 
mean  wind  speed  and  a  typical  direction.  Wind  stress  is  in  dynes/cm:;  heating  is  in  W7m:. 
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Figure  18.  Mixed  layer  deepening  is  calculated  with  time-dependent  forcing  from 
meteorological  data,  and  is  plotted  above  graphs  of  wind  speed  and  heat  flux.  The  periods 
of  more  rapid  deepening  correspond  to  the  periods  of  more  intense  forcing.  In  this  model  run 
the  heat  fluxes  were  corrected  by  scaling. 
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and  strongest  storm,  but  happened  gradually.  Even  though  the  final  forcing  event  was 
significantly  larger  than  the  other  two,  the  layer  was  more  resistant  to  deepening  because  it 
was  already  deeper.  For  a  deeper  layer,  more  energy  is  required  for  the  same  amount  of 
entrainment  simply  because  there  is  more  water  that  must  be  kept  turbulent. 

Although  the  magnitude  of  heat  removal  from  the  ocean  as  indicated  by  the 
meteorological  data  was  too  high,  there  is  a  high  degree  of  confidence  in  the  qualitative 
correctness  of  that  data;  the  storms"  occurances,  their  relative  magnitudes,  and  the  phase 
relationships  of  wind  peaks  to  heat  flux  maxima  are  taken  to  be  well  represented.  In  other 
words,  the  shape  of  the  surface  heat  flux  time  series  should  be  retained  even  though  its 
magnitude  should  be  adjusted.  Offsetting  the  heat  flux  by  uniformly  subtracting  52  W/m2 
results  in  the  correct  total  heat  loss  from  the  ocean,  but  results  in  periods  of  net  warming  of 
the  ocean  by  the  atmosphere  in  winter  which  are  thought  not  to  be  realistic.  Uniform  scaling 
(  multiplying  by    )  did  not  present  this  disadvantage,  and  so  was  preferred.  Another 

182 

option  of  first  adding  50  W/m2  as  an  approximate  additional  heat  loss  from  the  ocean  due  to 
net  long  wave  radiation,  then  multiplying  by  ,  was  investigated  later,  but  the 

182  +  50 

results  were  so  similar  as  to  be  almost  indistinguishable. 

Since  momentum  budgets  were  calculated  in  the  model,  they  provided  another  check 
on  the  reasonableness  of  one-dimensional  treatment  of  this  problem.  Integrating  the  mean 
wind-driven  U  and  V  of  the  mixed  layer  over  each  time  interval,  the  mean  wind-driven 
transport  is  computed.  The  progressive  vector  diagram  for  this  transport  is  shown  as  Figure 
(19),  and  indicates  the  water  column  would  move  ~5  km  to  the  west,  90  degrees  to  the  right 
of  the  mean  wind  direction,  which  is  almost  due  south. 
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Figure  19.  The  time  integral  of  the  wind-driven  velocities  U  and  V  of  the  mixed  layer  is 
plotted  as  a  progressive  vector  diagram,  showing  the  inertial  effects  and  the  mean  transport 
of  just  over  5  km  almost  due  west. 
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The  components  of  turbulent  kinetic  energy  axe  represented  by  u  - 2  ,  v 2  ,  and  w2, 
and  are  plotted  in  Figure  (20).  What  is  striking  in  that  figure  is  that  v  TKE  was  very- 
dominant  while  w  TKE  accounted  for  only  10  to  20%  of  the  total  TKE  during  storm  events: 
the  turbulence  was  very  anisotropic. 

Figure  (21)  shows  the  contributing  terms  to  u  TKE.  Shear  production  at  the  interface 
is  always  small,  and  pressure  redistribution  is  the  same  order  as  shear  production  at  the 
surface.  Figures  (22)  and  (23)  are  the  corresponding  plots  for  terms  contributing  to  the  v  and 
w  TKE  budgets.  In  Figure  (22),  shear  production  at  the  surface  is  large,  and  v  is  losing  even 
more  energy  through  pressure  redistribution  than  through  dissipation.  Figure  (23)  shows  that 
the  main  source  of  vertical  turbulent  energy  was  pressure  redistribution  of  wind  energy.  Also 
interesting  is  that  the  entrainment  buoyancy  flux  is  often  as  large  as  the  surface  buoyancy 
flux,  and  is  plainly  not  a  fixed  fraction  thereof*  The  integrated  effect  of  rotation  over  the 
period  is  neutral,  but  rotation's  effect  during  any  particular  storm  may  not  be.  Figure  (24) 
is  an  enlargement  of  one  section  of  the  time  series;  the  influence  of  rotation  during  the  first 
event  was  to  enhance  entrainment,  and  during  the  second  to  inhibit  entrainment.  These 
observations  show  that  for  the  same  magnitude  of  forcing,  a  different  depth  of  mixing  may 
result.  The  phase  relationship  between  wind  stirring  (via  pressure  redistribution)  and  surface 
buoyancy  flux  varies  for  different  events  as  well.  Could  this  also  be  a  factor  in  the 
effectiveness  of  a  particular  storm  for  deep  mixing? 


*  A  common  approximation  for  entrainment  by  free  convection  in  the  atmospheric  boundary  layer  is  that 
the  entrainment  bouyancy  flux  is  opposite  in  sign  and  a  fixed  fraction  (order  .2)  of  the  surface  buoyancy  flux 
(Stull,  1988). 
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Figure  20.  The  turbulent  kinetic  energy  components  are  plotted  as  functions  of  time  for  the 
time-dependent  forcing  case  with  surface  heat  flux  scaled.  Both  wind  stirring  and  cooling 
are  vigorous,  but  v  TKE  dominates  over  w  TKE. 
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Figure  21.  The  contributions  to  u  TKE  are  plotted  as  functions  of  time.  They  are:  shear 
production  at  the  surface,  shear  production  at  the  interface,  pressure  redistribution,  rotational 
redistribution,  dissipation,  and  storage. 
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Figure  22.  The  contributions  to  v  TKE  are  plotted  as  functions  of  time.  They  are:  shear 
production  at  the  surface,  shear  production  at  the  interface,  pressure  redistribution, 
dissipation,  and  storage. 
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Figure  23.  The  contributions  to  w  TKE  are  plotted  as  functions  of  time.  They  are:  surface 
buoyancy  flux,  entrainment  buoyancy  flux,  pressure  redistribution,  rotational  redistribution, 
dissipation,  and  storage. 
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Figure  24.  A  section  of  Figure  (23)  is  enlarged  to  show  detail. 
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Figure  (25)  shows  the  various  terms  contributing  to  the  total  TKE  of  the  mixed  layer. 
The  oscillatory  behavior  of  the  storage  term  is  an  interesting  depiction  of  the  mixed  layer's 
ability  to  alternately  store  and  then  dissipate  energy  in  response  to  wind  surges. 

Turning  to  the  question  of  why  the  final  mixed  layer  depth  on  1 9  March  is  -200  m 
too  deep,  two  basic  but  related  factors  are  felt  to  be  responsible:  the  model  does  not  retreat, 
and  it  may  transport  energy  to  the  bottom  of  the  layer  too  efficiently.  After  a  period  of  low 
winds,  when  the  wind  stress  suddenly  increases,  the  ocean  mixed  layer  requires  time  not  only 
to  re-stir  the  layer  before  entrainment  recommences,  but  to  transport  energy  surges 
downward  before  entrainment  velocity  can  increase.  Also,  just  after  increased  stirring, 
energy  that  is  concentrated  near  the  surface  of  the  ocean  gives  more  dissipation  than  the  same 
amount  of  energy  would  if  it  were  distributed  more  evenly  throughout  the  layer.  All  of  these 
effects  can  be  corrected  for  by  adjusting  the  model  constants:  the  dissipation  constant  m1?  the 
pressure  redistribution  constant  m2,  the  shear  production  constant  m3,  and  the  entrainment 
efficiency  constant  m4.  However,  in  this  research,  model  constants  were  not  adjusted  to 
achieve  quantitative  agreement  since  such  an  adjustment  would  have  meant  tuning  a  very 
general  model  to  match  one  specific  observed  case. 
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Figure  25.  The  contributions  to  total  TKE  are  plotted  as  functions  of  time.  They  are:  surface 
buoyancy  flux,  entrainment  buoyancy  flux,  surface  shear  production,  dissipation,  and  storage. 
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VII.  CONCLUSIONS 
A.         SUMMARY 

Coupling  of  the  ocean  to  the  atmosphere  through  deep  mixed  layers  in  polar  regions 
is  fundamental  to  our  climate  system,  and  much  of  the  energy  exchange  occurs  episodically 
during  storms.  In  order  to  understand  the  physics  of  the  ocean  mixed  layer  response  to 
storms,  a  bulk  turbulence-closure  mixed  layer  model  was  developed  and  was  used  to 
simulate  an  observed  case  of  deep  polar  sea  mixing. 

The  model  included  budgets  for  potential  temperature,  salinity,  horizontal 
momentum,  and  unsteady  three-component  turbulent  kinetic  energy.  In  addition  to  shear 
production,  surface  buoyancy  flux,  entrainment  buoyancy  flux,  pressure  redistribution,  and 
dissipation,  realistic  treatment  was  also  given  to  thermobaric  enhancement  of  buoyancy 
flux  and  to  Coriolis'  effects  on  turbulence.  The  model  demonstrated  the  ability  to 
qualitatively  simulate  events  when  initialized  with  actual  CTD  data  and  forced  with  a 
realistic  meteorological  time  series,  provided  the  heat  fluxes  were  adjusted  to  give  a  correct 
total  heat  loss  from  the  ocean.  Further  work  toward  tuning  of  the  model  constants  n\,  m2, 
m3,  and  m4,  and  the  inclusion  of  a  shallowing  mode  have  the  prospect  of  making  the  model 
quantitatively  useful. 

Several  conclusions  can  be  made  about  the  mixing  that  occurred  in  the  Greenland 
Sea  during  late  winter  1994: 
•  Horizontal  wind-driven  energy  exceeded  vertical  free  convective  energy  throughout 

the  period;  the  turbulence  was  anisotropic,  and  wind  stirring  was  dominant  over 

cooling  even  during  a  period  of  strong  cooling. 
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•  In  the  vertical  energy  budget,  pressure  redistribution  of  horizontal  energy  (from  the 
wind)  was  a  larger  source  of  energy  than  surface  buoyancy  flux;  free  convection 
was  not  dominant  even  when  considering  w2  only.  This  underscores  the 
importance  of  including  wind  unsteadiness. 

•  The  deepening  response  for  a  given  forcing  event  depends  not  only  upon  the 
strength  of  the  event,  but  also  upon  the  initial  turbulent  energy,  the  depth  of  the 
layer,  the  buoyancy  jump  at  the  bottom  of  the  layer,  and  the  stratification  below 
the  turbulent  region. 

•  The  direction  of  the  wind  may  be  important  not  only  during  strong  events,  but 
between  them  as  well,  because  it  influences  the  "spinup"  time  for  the  next 
entraining  event. 

•  The  phase  of  the  surface  buoyancy  flux  to  shear  production  also  differs  among 
events;  it  is  possible  that  this  difference  may  influence  the  storm's  efficiency  at 
deepening  the  mixed  layer. 

This  study  demonstrates  the  importance  of  water  column  initial  conditions,  the 
accurate  representation  of  strong  surface  cooling  events,  and  inclusion  of  the  thermobaric 
effect  on  buoyancy  to  determine  the  depth  of  mixing  and  ultimately  the  heat  and  salt  flux 
into  the  deep  ocean. 
B.         RECOMMENDATIONS 

One  of  the  primary  motivations  for  the  study  of  deep  mixed  layers  is  to  be 
eventually  able  to  parameterize  mixed  layer  dynamics,  in  particular  the  dynamics  of  deep 
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mixed  layers  and  subgrid  scale  penetrative  convection.  The  requirement  is  that  the  effects 
of  these  processes  be  included  in  general  circulation  models  despite  the  impracticability  of 
modeling  the  actual  physical  processes  within  the  general  circulation  model  itself.  This 
research  provides  a  tool  for  studying  various  forcing  events  and  their  calculated  effects, 
to  conceptualize  and  test  what  the  important  variables  in  such  a  parameterization  might  be. 
To  increase  this  model's  usefulness,  several  areas  need  further  study: 

•  Include  a  shallowing  mode  capability.  Turbulence  may  be  dissipated  even  without 
warming,  so  that  the  fully  turbulent  layer  may  at  times  be  much  shallower  than  the 
layer  that  is  well-mixed  in  potential  temperature  and  salinity.  The  additional  time 
required  to  re-stir  the  layer  is  currently  neglected. 

•  Include  a  finite  thickness  entrainment  zone.  Stability  within  or  across  the  zone  may 
be  important. 

•  Compare  deepening  and  energy  component  behavior  to  results  from  large  eddy 
simulations  of  the  same  events.  Specifically,  find  better  estimates  of  model 
constants,  and  ensure  the  correct  amount  of  anisotropy.  This  would  be  an  important 
step  to  improving  the  quantitative  accuracy  of  the  model. 
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